###############################################################
## Marginal Likelihood Computation
## Author: Xun Pang, Barry Friedman, Andrew Martin, Kevin Quinn
## Note: the function is to estimate the marginal likelihood
##       using the posterior means as the fixed their values of
##       the parameters in estimation
###############################################################


   Marginal.likelihood <- function(MCMCoutput,TimeIndex, beta0, B0, M, a00,
                                   b00, iteration, burninstage, y, X)
   {
    m <- nrow(MCMCoutput[[1]])
    beta.star <- apply(MCMCoutput[[1]],2, mean)
    if(M<=2) {
      bigP.star <- mean(MCMCoutput[[3]])
     }else{
       bigP.star <- apply(MCMCoutput[[3]], 2, mean)
     }
    state.mcmc <- MCMCoutput[[4]]
    S.start <- MCMCoutput[[4]][m,]
    ystar.mcmc <- MCMCoutput[[2]]
    ystar.start <- MCMCoutput[[2]][m,]

    # Prior ordinates
    Beta <- matrix(beta.star, nrow=M, byrow=TRUE)
    beta.prior <- sapply(c(1:M),
             function(k){dmvnorm(Beta[k,], beta0, B0, log=TRUE)})
    beta.prior.log <- sum(beta.prior)

    P.prior <- sapply(c(1:(M-1)),
             function(k){dbeta(bigP.star[k], a00,b00, log=TRUE)})
    P.prior.log <- sum(P.prior)
    cat("finish computing the prior ordinates, and reduced run begins","\n")
    
  ##################################
  ## Likelihood Ordinate
  ##################################
  # One reduced run needed here
   source("ChgPtProbitReducedMCMC.R")
   reduced <- Reduced.Run1(y=y,X=X,TimeIndex=TimeIndex, n=iteration,B0=B0,
                           burnin=burninstage,beta0=beta0, ystar.start=ystar.start,
                           beta.start=beta.star, S.start=S.start, M=M,
                           P.fixed=bigP.star,  likeihoodComputed=TRUE)
   cat("finish the reduced run, and compute the likelihood ordinate" ,"\n")

   likelihood.mean.log <- log(mean(exp(reduced[[1]])))
   cat("finish computing the likelihood ordinate, begin to estimate
        the posterior ordiates","\n")

   # posterior ordinates

    mcmc <-length(reduced[[1]])
    Beta <- matrix(beta.star, nrow=M, byrow=TRUE)
   
    beta.ordinate <- rep(NA, m)
    for (j in 1:m){
      state.current <- state.mcmc[j,]
      ystar.current <- ystar.mcmc[j,]
      betai <- rep(NA, M)
      for (i in 1:M){
         beta.state <- Beta[i,] 
         group <- which(state.current==i)

         TT <- length(group)
         if (TT==0) betadrawn <- rmvnorm(1, beta0, B0)
         else{
           XX <- data.frame(X)[TimeIndex%in%group,]
           ystar.M <- ystar.current[TimeIndex%in%group]
           NN <- length(ystar.M)
           total<-0
           for (k in 1:NN) {
             total<-total+outerself(as.matrix(XX[k,])[1,])
           }

           B1<-solve(total+solve(B0))
           beta1<-B1%*%(t(XX)%*%ystar.M+solve(B0)%*%beta0)
           betadprob <- dmvnorm(beta.state, beta1, B1)
         }
          betai[i] <-betadprob
        }
        beta.ordinate[j] <- prod(betai)
       }

     beta.mean <- mean(beta.ordinate)
     beta.ordinate.log <- log(beta.mean)

    # posterior ordinate of transition matrix
    reduced2 <- Reduced.Run1(y=y,X=X,TimeIndex=TimeIndex, n=iteration,B0=B0,
                             burnin=burninstage,beta0=beta0, ystar.start=ystar.start,
                             beta.start=beta.star, betafixed=FALSE, S.start=S.start,
                             M=M, P.fixed=bigP.star,likeihoodComputed=FALSE)
    state.reduced <- reduced2[[1]]
    pp.reduced <- c()
    for (q in 1:mcmc){
      state.current <- state.reduced[i,]
      stay.state <-  switchg(state.current, M=M)  
      bigP.prob <- c()
      for (i in 1:(M-1)){
       bigP.prob[i] <- dbeta(bigP.star[i], (a00+stay.state[i]),
                             (b00+1), log=TRUE)
      }
      pp.reduced[q] <- exp(sum(bigP.prob))
    }
    P.ordinate.log <- log(mean(pp.reduced))

    # Marginal Likelihood

    MargLik <- likelihood.mean.log+beta.prior.log+P.prior.log-
               beta.ordinate.log-P.ordinate.log
    cat("finish computing the marginal likelihood","\n")

    return(list(MarginalLikelihood=MargLik, LogLikelihood=likelihood.mean.log,
                Log.prior.beta=beta.prior.log, Log.prior.P=P.prior.log,
                Log.posterior.beta=beta.ordinate.log,
                Log.posterior.P=P.ordinate.log))
  }
